Analysis of the electric field gradient in the perovskites SrTiO,3 and BaTiOa: density 

functional and model calculations 
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We analyze recent measurements [R. Blinc, V. V. Laguta, B. Zalar, M. Itoh and H. Krakauer, J. 
Phys. : Condens. Matter 20, 085204 (2008)] of the electric field gradient on the oxygen site in the 
perovskites SrTiOs and BaTiOs, which revealed, in agreement with calculations, a large difference 
in the EFG for these two compounds. In order to analyze the origin of this difference, we have 
performed density functional electronic structure calculations within the local-orbital scheme FPLO. 
Our analysis yields the counter-intuitive behavior that the EFG increases upon lattice expansion. 
Applying the standard model for perovskites, the effective two-level p-d Hamiltonian, can not explain 
the observed behavior. In order to describe the EFG dependence correctly, a model beyond this 
usually sufficient p-d Hamiltonian is needed. We demonstrate that the counter-intuitive increase of 
the EFG upon lattice expansion can be explained by a s-p-d model, containing the contribution of 
the oxygen 2s states to the crystal field on the Ti site. The proposed model extension is of general 
relevance for all related transition metal oxides with similar crystal structure. 

PACS numbers: 77.84.DY, 76.60.-k, 77.80.-e 
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I. INTRODUCTION 

Perovskite compounds ABO3, with A being an alkali, 
alkaline earth or rare earth metal and B a transition 
metal element, attract much attention because of their 
importance both for fundamental science and for tech- 
nological applications^. Although the high-temperature 
cubic phase has a very simple crystal structure, this does 
not prevent these compounds to exhibit a large variety 
of physical properties rendering the perovskites to model 
compounds for studies of a large variety of different phys- 
ical phenomena. Within the perovskite family, we find 
superconductivity, e.g. in K^Bai-^BiO^, giant magne- 
toresistance, e.g. in LaMnO^, orbital ordering, e.g. in 
YTiO^ and ferroelectricity, e.g. in BaTiOa 7 -. The latter 
phenomena are of large interest because of technological 
applications. 

The compounds SrTi0 3 (STO) and BaTi0 3 (BTO) 
are usually considered to be isovalent. The valence and 
conduction bands of the two perovskites are formed by 
p-states of oxygen and d-states of titanium. In the high- 
temperature cubic phase, the Ti and O sub-lattices have 
the identical geometry for STO and BTO, the lattice 
parameters being a=3.8996 A£ and a=4.009 A-i respec- 
tively. As the temperature lowers, both compounds ex- 
perience a softening of an optical phonon mode, which 
corresponds to Ti motion towards the oxygen^. BTO 
exhibits a succession of phase transitions, from the high- 
temperature cubic perovskite phase to ferroelectric struc- 
tures with tetragonal, orthorhombic and rhombohedral 
symmetry^. In contrast, STO behaves as an incipient fer- 
roelectric in the sense that it remains paraelectric down 
to the lowest temperatures, exhibiting nevertheless a very 
large static dielectric response. It undergoes an antifer- 
rodistortive phase transition at 105 K to a tetragonal 



(I4/mcm) phase, but this transition is of non-polar char- 
acter and has little influence on the dielectric properties^. 

The first determination of the 17 O electric field gradi- 
ent (EFG) on the oxygen site in perovskites was recently 
reported for STO and BTO^ together with first-principle 
calculations using the linearized augmented plane wave 
(LAPW) method was used. The most striking feature 
in the experimental and theoretical data is the large dif- 
ference of the EFGs between the two compounds. The 
calculational investigation of Blinc et al. concluded, that 
the mag nitude of the EFG of 17 O in BTO is larger than 
the EFG of 17 in STO due to two effects: (i) larger 
lattice parameters in BTO compared to STO and (it) 
a larger ionic radius of Ba compared to Sr. While the 
experimental determination (NMR) can not provide the 
sign of the EFG, the LAPW calculation yielded a nega- 
tive EFG. A negative EFG corresponds to a prolate elec- 
tron density, which implies the importance of covalence 
effects. 

In order to elucidate the origin of the sign of and 
the different contributions to the EFG, we have per- 
formed first-principle calculations using a local orbital 
code (FPLO^) that is especially suited to address these 
questions due to its representation of the potential and 
the density allowing easy decomposition. The calcula- 
tional details of our investigation are given in Sec. [HI 
and the obtained results are presented in Sec. IIIII These 
results can not be explained by intuitive models, which 
are also described in this section. Therefore, a more 
complex microscopic model Hamiltonian is introduced in 
Sec. IIV1 Using the properties of this p-d like Hamilto- 
nian, an agreement with the obtained experimental and 
theoretical results and a deeper, microscopically based 
understanding is obtained. 
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Table I: The experimental and calculated values of the EFG 
(in 10 21 V/m 2 ) on the oxygen site in the cubic phase of the 
two perovskites. The last four lines refer to equations given 
in the appendix. 



II. CALCULATION METHODS 

The electronic band structure calculations were per- 
formed with the full-potential local-orbital minimum 
basis code FPLO (version 5.00-19)^ within the lo- 
cal density approximation. In the scalar relativistic 
calculations the exchange and correlation potential of 
Perdew and Wangi^ was employed. As basis sets 
Ba (4d5s5p/6s6p5d+4f7s7p), Sr (4s4p/5s5p4d+6s6p), Ti 
(3s3p3d/4s4p4d+5s5p) and O (2s2p3d+3s3p) were cho- 
sen for semicore/valence+polarization states. The high 
lying states improve the completeness of the basis which 
is especially important for accurate EFG calculations. 
The lower lying states were treated fully relativistic as 
core states. A well converged A:-mesh of 455 fc-points was 
used in the irreducible part of the Brillouin zone. 



III. FPLO ANALYSIS RESULTS 

In FPLO, the EFG on a nucleus at a given lattice site 
may be represented as the sum of two contributions: An 
on-site contribution V™ (see Eq. (|27|) ). which comes from 
the on-site contribution of the electron density of the 
given lattice site, and a second term, the off-site contribu- 
tion V°P (see Eq. 1(28)1 ). which results from the potential 
of all other atoms (see App. A). The on-site contribution 
V™ can be analyzed further. It can be split up in p-p, 
s-d and d-d contributions (see App. B). 

The on- and off-site contributions as well as their sum 
and the dominating p-p contribution (see Eq. ((3TJ) ) are 
shown in Tab. [B Whereas the total EFG for 17 in 
BTO agrees well with the experiment (1 % deviation), 
the total EFG for 17 in STO is in discrepancy with 
the experiment (38 % deviation), see Tab.|U Compared 
to the EFGs calculated with the LAPW code in Ref. 0, 
we obtain almost the same absolute value of V zz but the 
opposite sign, see Tab. HI Our calculated EFGs as a func- 
tion of the lattice parameter a for both compounds reveal 
the same tendency as observed in Ref. 0: The absolute 
value of the EFG increases under the lattice expansion 
(see Fig. [J). From Fig. Q] we also conclude that the EFG 
of BTO is not only larger than the EFG of STO due to 



larger lattice parameters ("lattice effect") , but also due to 
an "cation effect", which is responsible for the remaining 
difference. This lattice effect is demonstrated by the shift 
between the two EFG curves in Fig. [TJ 




lattice parameter a [A] 



Figure 1: Calculated V zz in dependence of the lattice parame- 
ter a. Vzz for the experimental lattice parameter is marked by 
a shaded square. The "cation" and "lattice effect", which are 
responsible for the difference in Vzz for these two compounds 
are indicated by the red and black arrow, respectively. In- 
set: The anisotropy count Ap (see text) in dependence of the 
lattice parameter a. 

The increase of the (absolute value of the) EFG upon 
lattice expansion is rather counter-intuitive. In the tra- 
ditional approach, the spherically symmetric electronic 
shell of an ion is perturbed by the potential of the ex- 
ternal (point) charges of the solid. As a result, the total 
EFG on the ion nucleus is caused by the EFG of the ex- 
ternal potential, and is roughly proportional to it. It is 
clear that this approach predicts the opposite tendency: 
The strength of the external potential is inversely propor- 
tional to the lattice constant and thus the (absolute value 
of the) EFG should diminish under the lattice expansion. 
The failure of this approach to describe the observed be- 
havior of the EFG indicates that a fully ionic description 
of the perovskites is inappropriate. 

In an alternative approach, the electronic shell of the 
atom is disturbed by the hybridization of the wave func- 
tions with the states of the surrounding atoms. The 
hybridization results in the asymmetry of the electronic 
cloud of the atom and the EFG on its nucleus. Appar- 
ently, this covalent approach predicts the same tendency 
as the ionic one: It is usually believed that the hybridiza- 
tion diminishes with the increase of the bond length. In 
both approaches we may say: When expanding the lat- 
tice, we diminish its influence on the atom, and the elec- 
tronic shell should become closer to that of the free atom. 
Hence, we come to the conclusion: The (absolute value 
of the) EFG should diminish under the lattice expansion, 
which is opposite to the experimental observation and the 
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Figure 2: The on-site V™, off-site V°/ f and total EFG as a 
function of the lattice parameter a. The grey shaded squares 
mark the experimental lattice parameter for V zz . 



results of both first-principle calculations. We will tackle 
this problem in detail in Sec. HVl 

Another problem it the different sign of the EFG ob- 
tained from the two different band structure codes. If 
the sign of the EFG is taken into account, the slope in 
our graph (Fig. [I]) is opposite to the slope in the graph 
obtained with the LAPW code (Fig. 5 in Ref. H). Since 
the NMR experiment is not sensitive to the sign of the 
EFG, we will investigate the influence of the lattice ex- 
pansion on the different contributions to the EFG to get 
more insight in this issue. 

Our calculations show that both the on-site and the 
off-site contribution to the EFG have comparable val- 
ues for the perovskite lattice, see Tab. U and Fig. O In 
Fig. [21 the two contributions, V™ (dashed line) and V°J f 
(dash-point line) and the total EFG (full line) are shown. 
Whereas the off-site EFG decreases only slightly upon 
lattice expansion, the on-site EFG increases strongly with 
increasing lattice parameters, resulting in the significant 
increase of the total EFG. We also observe that the off- 
site EFG is almost identical for these two structures, 
which is in line with the observed weak dependence of 
V°P on the lattice parameters. 

The on-site EFG is mainly caused by electrons with 
p character, see table |H Therefore, we will investigate 
the corresponding anisotropy count Ap^. In the per- 
ovskite structure AB0 3 , the oxygen site has axial sym- 
metry, and the z-axis is directed along the B-0 bond. 
Thus, the anisotropy count is the difference between the 
population of the oxygen 2p a- (corresponds to p z ) and 
7T- (corresponds to p x , y ) orbitals. In the inset of Fig. [T] we 
see that the anisotropy count Ap increases with the lat- 
tice expansion. This is in agreement with the increasing 
on-site EFG. If we focus on BTO, where the experimen- 
tal and calculated (for the experimental lattice parameter 
a = 4.009 A) value for the EFG agree very well, we see 



Figure 3: Occupation of p x and p z states in dependence of 
the lattice parameter a. 



that this positive V zz corresponds to a positive Ap. That 
means the p electron density (responsible for the EFG) 
has an oblate shape, since more electrons are occupying 
the p^-orbitals than the p z -orbital, which is in agree- 
ment with the positive sign of the EFG. 

After concluding that the sign of V zz for 17 for both 
STO and BTO should be positive, we come back to the 
counter intuitive behavior of the increasing EFG upon 
lattice expansion. Fig. [3] reveals that the increase of Ap 
under lattice expansion, which is responsible for the in- 
creasing EFG upon lattice expansion, is due to an in- 
creasing occupation of pi- (corresponds to p x ,y) and an 
decreasing population of a- (corresponds to p z ) orbitals. 



IV. DISCUSSION 

In order to understand this anomalous behavior of 
the cr-orbital, we will analyze the main features of the 
electronic structure of perovskites. Detailed band struc- 
ture studies of perovskite compounds were performed by 
Mattheis a 12 i 13 i 14 , who also proposed a first tight-binding 
fit for the band dispersions. Wolfram et a / , 15 i 16 i 17 ( c f. 
also Ref. [T3) developed a very simple model (Wolfram 
and Ellialtioglu, WE) for the valence and conduction 
bands, which reflects their basic properties. The WE 
model includes the d-orbitals of the B ion and the p- 
orbitals of the oxygen. Wolfram et al. pointed a quasi- 
two-dimensional character of the bands out, which is due 
to the symmetry of the orbitals. If one retains only 
nearest neighbor hoppings, the total 14 x 14 Hamilto- 
nian matrix (five ri-orbitals and 9 p-orbitals) acquires 
block-diagonal form at every value of the momentum. 
The three 3x3 matrices describe the 7Ty -bands (ij = 
xy,yz,xz). Every dy -orbital of the t 2g symmetry cou- 
ples with its own combination of oxygen 2p 7r-orbitals, 
which lie in the same plane perpendicular to the bond 



direction. They form a pair of bonding and anti-bonding 
states. The remaining combination of the 2p 7r-orbitals 
in the same plane forms the non-bonding band. Wol- 
fram et al. call this group of bands 7r-bands. The states 
described by the 5x5 block matrix are called cr-bands, 
since they are formed by oxygen 2p er-orbitals, which are 
coupled with the e g (d x 2_ y 2 and d z 2) orbitals of the B 
ion. This matrix decouples into one non-bonding band 
and two pairs of bonding and anti-bonding bands. 

Fig. [4] shows the calculated band structure for STO 
for two different lattice parameters a. The features men- 
tioned above are clearly seen (cf. Fig. 2 of Ref. U3). The 
anti-bonding 7ry-bands are situated between 2 and 4 eV, 
where the 7rj, 2 -band is almost dispersionless in the direc- 
tion r — > X. This manifests the quasi-two-dimensional 
character of the bands. The bands originating from the 
d e g -orbitals are in the range from 4 to 8 eV, where the 
band expressing d z 2 character is dispersionless along the 
r — > X direction. The valence band has a more com- 
plex character due to additional mixing from the direct 
p — p hopping. This is neglected in the simple version 
of the WE model. Nevertheless, we see that the non- 
bonding bands lie on top of the valence band and have a 
much smaller dispersion than the bonding bands, which 
lie below — 1 eV (^ij) and below — 3 eV (cr-bands). The 
latter have a larger dispersion due to much larger d — p 
hoppings. 

Although the Kohn Sham theory is not good for ex- 
citation spectra, or obtaining the correct energy gap, it 
yields reliable occupation numbers, on-site energies and 
transfer integrals, especially in the absence of strong cor- 
relations. Therefore, we can use our LDA band structure 
to obtain reliable parameters as input for further treat- 
ment using model Hamiltonians. 

In the following, we explore within the WE model 
how the occupation numbers and the resulting anisotropy 
count for the p-orbitals depend on the lattice parameters. 
In dielectric compounds like STO and BTO, the bond- 
ing and non-bonding states are fully occupied. Contrary 
to the non-bonding bands, which have almost pure p- 
character, the bonding and anti-bonding bands are mixed 
p-d-bands. The population of the p-orbitals is given by 
the sum of the occupation numbers of the non-bonding 
and the bonding bands, whereof the latter are lattice pa- 
rameter dependent. 

Every pair of bonding and anti-bonding states is de- 
scribed by an effective two-level modeU& 

H m = A m (|d,k)(d,k|-|p,k)( P ,k|) 

+V m / mk (|d,k><p,k| + |p,k><d,k|). (1) 

Here, m describes the character of the band m = tt, a 
and /mk is a dimensionless function, which depends on 
the dimensionless variable ka (note that k is measured 
in units of n/a, so neither ka nor / m k depends on a). 
The state mixing is defined by the interplay of the on- 
site energy difference A m and the transfer integral V m , 
which determines the bandwidth of the corresponding 
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Figure 4: SrTiO;?: band structure for two different lat- 
tice parameters a = 3.8996 A (black/colored full lines) and 
a = 4.009 A (brown dashed lines). The different band char- 
acters are given by different colors: blue (bonding, cr), cyan 
(bonding, n), orange (non-bonding), red (anti-bonding, mj) 
and green (anti-bonding, d eg ), see text. Since it is not easy 
to interpret the valence band, the colors in the valence band 
are only approximate. 

band. The eigenstates of the Hamiltonian Eq. ([I]) have 
the form 

|k, v) = Cd-kv d, k) + c pkl/ \p, k) , (2) 

and for the p-states, the following energies and occupa- 
tion numbers are obtained 

E km u = iy^Al + V r l(f m]i ) 2 , v = ±l (3) 
npkmv = 2 Icpk^l 2 = 1 — A m /Ek v . (4) 
ndkmv = 2 \cdkv\ 2 = 2 (l - |c p k^| 2 ) = 1 + -pr- (5) 

Here, v = +1 describes the anti-bonding and v = — 1 de- 
scribes the bonding band. In this two-level system, two 
asymptotic behaviors are possible. First, A m /V m — > oo, 
which yields for the occupation numbers of the bonding 
bands n p km,-i ~* 2 and ndkm,-x — * 0. In this case, both 
electrons are in the p-state of the ligand ion and d-states 
are empty, called the ionic limit. Second, A m /V m — > 0, 
which yields for the occupation numbers n p km,-i - * 1 
and nscm-i — * 1- In this case, the electrons are equally 
shared by the p-, and d-states. This is the covalent limit. 
From the trends in Fig. [3l we observe that while the pop- 
ulation of the pTr-orbitals increases, the population of the 
Pa-orbitals decreases. This means the Ti-0 7r-bond gets 
more ionic under lattice expansion (as expected) whereas 
the Ti-0 cr-bond gets more covalent, which we will try to 
explain with this model. 

The parameters of this model can be extracted from 
the band energies at symmetry points of the Brillouin 
zone in Fig. [4] (see the App. C for more details). 
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For example, the on-site energies A m can be obtained 
from the L point, since due to symmetry, the d—p mixing 
vanishes at this point and the band states acquire a pure 
d or p character. For a = 3.8996 A we have for STO 
Ed t2g « 1-7 eV, E dcg w 4.3 eV, E p ps -1.2 eV. This yields 
(using Eq. {42J) and Eq. gTJ) 2A T = £ df2s -£! p w 2.9 eV 
and 2A ff = E dsg - E p ps 5.5 eV. 

From these values and the / m fc as given in Refs. 
we obtain the Slater-Koster hopping parameters V a ps 
2.1 eV Eq. {44]), K = V p(i7r ps 1.6 eV Eq. {43]) and « 
2.7 eV Eq. ijgfijl^. 

Since the occupation numbers of the non-bonding 
bands do not depend on the lattice and the anti-bonding 
bands (u = +1) are not occupied, we consider the bond- 
ing bands (v = —1) only. The contributions from the 
bonding bands to the population of the p m orbitals n Pm 
are obtained by a sum over the Brillouin zone. In order 
to analyze the occupation in dependence of the lattice ex- 
pansion, we need the derivative of the occupation number 
with respect to the lattice parameter a. From Eq. (J4j) we 
obtain for the derivative (denoted by /) 



^(/ mk ) 2 A, 



A? n + vi (f mk y 



A' 

m 

A m 



V , 

f- ) • (6) 

1 rn 



The derivative of n Vrn is proportional to 



A m 



V 

v m 

v m 



(7) 



Fig. [3] shows that n' Pm has a different behavior for m = a 
(n' Pa is negative) and m = ir (n' p ^ is positive). Thus, 
within the WE model, the observed increase of the EFG, 
which is due to the decreasing occupation of the p a - 
orbitals would yield 



and strongly a-dependent contribution for m = a from 
the semi-core s-states of the ligand. Indeed, (cf. Fig. |4]) 
the change due to the increasing lattice parameter a is 
much larger for A^ than for A^. The main electrostatic 
contribution, which implies ScF.ei cx a , leads to 



A' S, 



CF.el 



Since e d - e p + S C F,ei > 5 C F,ei is S C F,ei/^a < 1 and 
therefore 



A,, 



(10) 



Combining the estimates from Eq. (J9j) and Eq. (fT0)| . we 

get 



A' 



< 1 < 3 < 



V 

v a 



(ii) 



This is in contradiction to the inequality ([8]), leading to 
the conclusion that the WE model, though consistent 
with the intuitive expectations (see Sec. IIII|1 is unable to 
predict the observed behavior of the tr-orbital occupation 
in Fig. 03 

A possible reason for the failure of the WE model is 
that according to Ref. a large contribution to the CF 
comes from the oxygen 2s-orbitals, which lie almost 18 eV 
below the Ti 3d level, A sd = 17.9 eV, but have a large 
matrix element V sda = 3.0 eV with the e g orbitals. This 
suggests to extend the WE model taking into account the 
oxygen 2s-states in order to explain the increasing EFG 
upon lattice expansion. This is Harrison's model, where 
V sda is obtained from Eq. (|45|) 



r i2 



£s + £d ± x l(- s 



e d 



) 2 + gv^ 



)1 

V a 



< 



Act 



(8) 



Both Ao- and V a decrease upon lattice expansion: 
Fig. |4] shows that the energies at the T-point E d<ig and 
E dt2g and the bandwidths are smaller for the larger lattice 
parameter a = 4.009 A, than for the smaller lattice pa- 
rameter a = 3.8996 A. A commonly accepted estimate^ 
for the dependence of hopping integrals on a is Va oc a~ a 
with a between 3.5 and 4 (from the LDA band structure, 
we obtain a = 3.5 ± 0.5). This gives 



V' 

a-j- — a > 3. 

v a 



(9) 



Act is the difference in energy of the atomic levels cor- 
rected by the crystal field (CF)2& A a — s d ~ e p + ScF,a- 
The crystal field consists of two contributions^: A (dom- 
inating) electrostatic contribution, which is the differ- 
ence of the Madelung potentials of Ti and O, hence 
ScF,ei a_1 , and a hybridization contribution, which, 
in our case (octahedral coordination), contains a large 



with e s = —16.2 eV, e d = 1.7 eV and Ti 2 = 4.3 taken 
from the band structure. 

Taking the s-orbitals into account, V a in the inequal- 
ity JSJ is replaced by V sda . Harrison^ argues that the 
a dependence of V sda is similar to the a dependence of 
V pd a- This suggestion is confirmed by our LDA calcula- 
tions. Thus, we obtain 



Vsda' 
V s da 



a 
(I 



(12) 



On the right hand side we have the on-site energy dif- 
ference, which is given by A CT ~ A^ + 3V^ d(J / A sd , cf. 
Eq. lj49j) . The derivative of this expression is 



A^ ~ -7 VsdaV' sd(J . 

<±sd 



(13) 



Note, that here we assumed A'^ — A' sd — 0. Applying 



Eqs. (JUl) and (O yields 



6Vi 



a A sd A w + 3V/ do 



(14) 



Inserting Eq. lfT2|) and Eq. lfT4|) in the inequality ([8]), we 
obtain within the Harrison model the observed increase 
of the EFG, due to the decreasing occupation of the p a - 
orbitals, if the following inequality is fulfilled: 

« = _W l K _a 6V 2 da 

a V sda A a a A sd A v + W 2 da 

^A sd A n < V 2 da . (15) 



lattice expansion, see Fig. O Also the charge redistri- 
bution may change the value of V°/' , but as we see in 
Fig. El it has a minor effect: The off-site EFG for BTO is 
smaller than for STO, but the distance between the two 
curves is smaller than the lattice parameter dependence 
of the two curves. 



Using the values obtained from the LDA band struc- 
ture {V sda = 3.0 eV, A sd = 17.9 eV and A T = 1.4), we 
see that Eq. lfl5|) is fulfilled. Thus, the inequality Eq. ^ 
holds for the STO er-orbitals and the observed negative 
slope of n z in Fig. [3] can be understood. 

After revealing the origin of the counter-intuitive be- 
havior of the on-site EFG, we will discuss the unusually 
large value of the off-site EFG of the considered com- 
pounds. The dependence of this contribution with re- 
spect to the lattice parameter can be estimated in the fol- 
lowing way: From the multipole expansion of a potential 
of a given ion, the sum of the monopole contributions to 
v°f*(r) Eq. f22|) has the slowest convergence. This con- 
tribution may be calculated within a point charge model 
(PCM). Therefore, we note that the V zz value created in 
the origin by a unit charge situated at the point R equals 
the value of the z-component of the electric field E z , cre- 
ated in the origin by the unit dipole directed along z-axis 
and situated at the same point R: V zz = (3Z 2 — R 2 )/R 5 . 
That means, for the calculation of the EFG within the 
PCM, we need the electric field S"(r) of dipoles located 
at the sites R, which are polarized along the z direc- 
tion and whose polarization is unity, at various points r 
through the cubic lattice: Six) = X)r-^(-^ — r )- Here, 
r = a(x, y, z) and R = a(l, m, n) with a being the lat- 
tice parameter and i, m, n = 0, ±1, ±2. Using Eq. (16) of 
Ref . l2ll we obtain for the EFG in the PCM at the oxygen 
site 



V, 



PCM 



+2n o 5(0,ii) 

= — % [30.080n Ti - 8.668 {n A - no)] ■ (16) 
a 6 

Here, riTi is the monopole moment of the ionicity of Ti. If 
we insert the charges of the Ti ion titi, the O ion no, and 
the A ion ua = — (n-n + 3no) (with A=Sr, Ba) obtained 
from the FPLO calculations, we obtain e.g. for STO 
Vzz GM = 1-30 • 10 21 V/m 2 . This value is very close to 
V° z ff = 1.19-10 21 V/m 2 , see Tab.H So, we obtain a good 
agreement for the EFGs obtained from the simple PCM 
model and the more complex calculation. This means, 
the FPLO code yields realistic relations of the charge 
distributions. 

The prefactor e/a 3 in Eq. fl6|) is responsible for the 
observed decrease of the off-site contribution in case of 



SUMMARY AND CONCLUSION 



In summary, we have performed first principle calcu- 
lations of the electric field gradient on the oxygen site 
for BaTi03 and SrTi03 for different lattice parameters 
a. The values of our calculated EFGs agree well with the 
measured and, apart from the sign, with the calculated 
(LAPW) counterparts from Ref. y|. 

Decomposition of the EFG yields a large on-site contri- 
bution originating from the oxygen 2p shell. The on-site 
EFG reveals an anomalous dependence of the Po-orbital 
population with respect to the lattice parameter a: The 
population decreases under lattice expansion, i.e. the p-d 
hybridization grows with increasing Ti-0 distance. Sim- 
ple ionic and covalent approaches lead to the conclusion 
that this behavior is counter-intuitive. Also the effective 
two-level Hamiltonian proposed by Wolfram and Ellial- 
tioglu, which describes the relevant states of the valence 
region (oxygen p- and titanium d-states) fails to describe 
the observed behavior of the EFG upon lattice expansion. 
Only the inclusion of the 2s states to the crystal field 
results in a consistent picture: In fact, lattice expansion 
causes a charge transfer from the p a - to the s-orbitals of 
oxygen, whereas the population of the oxygen 7r-orbitals 
increases with a. This charge redistribution leads to the 
increase of the EFG, which is the main reason for the sur- 
prisingly large difference of the EFGs between BaTi03 
and SrTi0 3 . 

We expect that the observed feature, the increase of 
the anisotropy count of the p-shell with the bond length, 
is common to all d-metal-oxygen bonds and should be 
taken into account accordingly in the interpretation of 
the relevant experiments. 

The considered ^4Ti03 systems are not strongly corre- 
lated, since the Ti 3c? shell is formally empty. For mag- 
netic ions with partially filled c?-shells, the influence of 
the 2s orbitals will be diminished because the charge 
transfer energy A sd will include the on-site Coulomb re- 
pulsion within the d-shell. 

As a side effect, our investigation sounds a a note of 
caution: When performing a mapping of a complex DFT 
band structure calculation onto a microscopically based 
minimal model in order to gain deeper physical under- 
standing, care has to be taken that all relevant interac- 
tions are included. 
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VI. APPENDIX 
A. EFG implementation in FPLO 

The EFG is a local property. It is a traceless sym- 
metric tensor of rank two, defined as the second partial 
derivative of the potential v(r) evaluated at the position 
of the nucleus 



d 2 v(r) 



1 



r=0 



With the definition 



lim — D 2m (r), 



(17) 



(18) 



uses the Ewald method to handle the long-range inter- 
actions (see^ section D). The density is modified with a 
Gaussian auxiliary density ni(r) = ni(r) — nf w (r)21. In- 



serting this modified density in the potentials Eq. d 
and Eq. ((22]) yields 



v(r) 



l (r) 



,Ew,on 



(r) 



6 o// (r) 



v 



Ew,ojf ( 



•)• (23) 



These contributions are calculated to get the total EFG. 

The first contribution is £>°"(r) in Eq. {231). This 
potential is given by Eq. IpTj) using the modified den- 
sity n SQ _L{r'). The corresponding v SOl 2m(r) components 
needed in Eq. (fl8|) are obtained from the solution of the 
radial Poisson equation (see Ref. Eq. (49)) 



v So ,L(r) 



47T 

21 + 1 



1 



dxx l+2 n So . L (x) 



+r l I dxx l+1 h So , L (x) 



Using the rule of L 'Hospital we obtain for the V^m com- 
ponent (from which is obtained) 



''so ,2m 



(0) 



dxx 1 n S0 2m{x) 



(24) 



can the Cartesian EFG tensor Eq. (fT7|) also be expressed 
in (real) spherical components (I = 2, m = ±2, ±1,0) 



V21 \ 



V 



V21 



■Vai 



Vo 



2,-1 



(19) 



In FPLO, the EFG on a nucleus at a given lattice site So 
may be represented as the sum of two contributions 



V i4 



v on {r) 
v off {r) 



, 3 ,n S0)L (\r'\)Y L (v') 



(21) 



^--I^-a) [v™(r)+v°ff(r)] (20) 

E / ** 

L J 

E 

E u 



d3r ,n s MW\)Y L (ri) 

r-R-s-r' v ; 



R+s^so 



where Yl are the (real) spherical harmonics; R is a Bra- 
vais vector, and s is an atom position in the unit cell. 
The index L = nlm also absorbs the spin and the prin- 
cipal quantum number. The first term in Eq. (|20]) . the 
on-site contribution, comes from the on-site contribution 
of the electron density of the site so, and the second term, 
the off-site contribution, comes from the potential of all 
other atoms. 

Since the angular momentum components of the local 
charge density give rise to multipole moments, which de- 
termine the Coulomb potential for large distances, FPLO 



The first term in Eq. lf24|) is the 2m component of the 
electronic density at the nucleus ro So ,2m(0) = fi SOi 2 m (0). 
The ri2m component of a spherical harmonic expansion 
of an analytic function around a given point behaves as 
n-2m — 0(r 2 ). The only non-analyticities of the electron 
density are caused by the spherical singularities of the 
nuclear potential and this can not be aspherical. There- 
fore 7i2 m (0) = 0, which can be shown explicitly both in 
a non-relativistic and full relativistic theory. 

The second contribution is v°^(r) in Eq. |23)l . This 
potential is given by Eq. (|22|) using the modified density 
n s .L{r'). Since the density n Sj 2m is not given at the site 
So, where the atom under consideration is sitting, this 
equation has to be expanded. This can be done explicitly 
but the derivation as well as the result for are very 
bulky 2 ^ and therefore not given here. 

The third contribution are v Ew >° n ( r ) + v Ew ' off (r) in 
Eq. (|23]). which have to be calculated from the Ewald 



density alone. The auxiliary density nf w (r) is given as 
a Fourier expansion, resulting in the Ewald potential in 



Fourier space 



Eu: 



„Ew 



, Eq. (52) in*. V z 



Ew 



IS 



obtained by differentiating v (r) = J2c 



3 iGs 1 ,Ew 



V t = - E W^) K(e* Gs vg™) (25) 



G 



The total EFG tensor Vij is given by the sum of these 
three contributions 



V- ■ — V on + V off -I- 



(26) 



In order to analyze the on-site and off-site contributions, 
we define the on-site EFG as being the first term in 



Eq. l(20|) . but calculated from the unmodified density 



dxx 1 n SOt2m (x). 



5 jo 

The off-site EFG then is taken to be 



V 2m — "2m ~ K 2m . 



(27) 



(28) 



a [A] 


r M 


Ti 


Tl5 


r 25 


r 15 


r 25 


r u 


X 5 


*i 


3.8996 


-17.199 


-16.177 


-2.891 


-1.166 


-0.372 


1.709 


4.319 


3.705 


6.551 


4.009 


-16.923 


-15.968 


-2.828 


-1.046 


-0.408 


1.579 


3.800 


3.332 


5.798 



Table II: The energies at the F and X points in SrTi03 given 



in eV. Here, Ti 

Ed. a 



F25 



T'25 = Ed, : 



Ed and Fi2 



I ! . Orbital contributions to the EFG 



C. Background for Sec. HVl 



In FPLO the electron density is separated into a net 
density and an overlap density (see Ref.0 section B). The 
dominating net density is calculated from two orbitals at 
the same site R + s = R' + s' = s 



net/ \ \ ^ k,n / \ *k,n / \ 

n *o ( r ) = c So L 1 ( Ps ,L 1 (r-s )-c so - L2 ip s ^ L2 (r-s ). 

k,nLi,L2 

The basis functions <p So ,L are localized on the lattice sites 



^so,i(r - s ) 



|r - s \)Y L (r - s ) 



The 2m component of the radial net density, needed for 
the contributions of the net EFG, can be calculated from 



% m {r) = J < e *(r)r 2m (r-s )^ 



(29) 

= E c ^^ (|r - So|)^(|r - so\)GTXr\ 

L1.L2 

where G™ 1 ^™ 2 '" 7, are the Gaunt coefficients and cl x Li = 

n c s Li c soL™ • D ue *° ^ ne properties of the Gaunt co- 
efficients, Tij o e ' 2m consists only of p-p, d-d, and s-d (and if 
present p-f and /-/) contributions. These contributions 
to the on-site net EFG V zz ,net are obtained by inserting 
Eq. 429} into Eq. E.g. the p-p contribution V^™;™^* 

is calculated from 

VZ'TJ = 2^ dxx^nZfZix) (30) 



2m, pp 

C&) = [^>)] 2 E C M m2G M 



ni2 ,m 
2 



The main component V°™j™ et = "j^%a,pp ^ s calculated 
from 



U^ (-)] 2 Y,( c t^i:o (3i) 

k,n V 



k,n *k,n . k,n *k,n 

's ,l,-l C s ,l,-l + C S ,1,1 C S ,1,1 , 



We see, that this density is proportional to the difference 
of occupation in p z (m = 0) and p XiV (m = ±1) states, 
which is the anisotropy count. 



In order to extract the parameters from the band struc- 
ture we need the total Hamiltonian 

^ = E [ffm + e m (|d,k)(d,k|4>,k)(p,k|)] . (32) 



Here, H m is the Hamiltonian given in Eq. ([I]) and e m is 
the mean energy of a pair of bands. The energies are 
therefore obtained from 



E km „ + e m = e m + u^] A 2 m + V^{f mk ) 2 . (33) 

For the three pairs of the 7Ty bands, / m k is given by 

fl jk = 2(2-Q-C 3 ) , (34) 

with C, = cos(fcia). The two tr-bands are distinguished 
by the index A = ±1 and / m k is 



/a A k — 3 C x C y C z 



+X {Cl + C 2 + C 2 Z - C x C y - C X C Z - C y C z ) 1/2 (35) 

Inserting these in Eq. I|33p for the T point (ka = 0), and 
the X point, (k x a = ir, k v = k z = 0) we obtain 



r 25 

and T 2 5 

r' 

1 25 

and r 25 

Tl2 

and Ti2 

X 5 



^7T &(T ^CJ; 

Ep ~ Sp 

E d t 2 !l ~ £ d 



(36) 
(37) 
(38) 

(39) 
(40) 



Here, e denotes the energy of the atomic level, and E 
denotes the energy level corrected by a 'crystal field' 5c f, 
see below. 

Now it is trivial to find the parameters A m , V m 



2 An — r 25 - T 2 5 = E dt2g - E p , 

2A a = T12 — T 25 = E c i cg — E p = 

4V 2 = (X 5 -T' 25 )(X 5 -T 25 ), 

w 2 = {x x - r 12 ) {x ± - r 25 ) . 



£d — £p 



(41) 

■Scf,(42) 
(43) 
(44) 
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The energy values at the different T and X points for 
SrTi0 3 are given Tab. HI 

So far, we have used the WE model, i.e. we have taken 
into account only the oxygen p and the titanium d states. 
Since this model is not sufficient to explain the observed 
behaviour of the oxygen Po-states, we have to expand 
the model. Harrison's modetf^ includes also the oxygen 
s-states. The s-states change the dispersion in the a 
bands, so that we have two parameters V p da, V s da in- 
stead of just one V a . Thus, the expressions become more 
complex, even in the symmetry points. In this model, 
the Eqs. l[36|) . (|37)l and (|39j) remain the same and the pa- 
rameters e p , Sd and are unchanged. For Ti2 Eq. lj38|) 
and Xi Eq. (J40J) , Harrison obtains 



Hence, 



Ed ca = Ti2 ~ Ed H — r~ ~ j 



For the main text, we need an expression for A CT : 

K = (ria-r 26 )/2 
= (r 12 - £p )/2 



1 



Ed 



6K 



2 V"" ' A sd ° p 
A v + 3V s 2 d JA sd . 



(49) 




(45) 



W p 2 dcr , (46) 



where Eda = £d + 2V^ CT /A s( j. From these equations, 
the parameters V pd a and V a da can be obtained. Besides, 
there is also an additional equation for Ti 



(47) 



Substituting A sd = Ed - £ s » V s do in Eq. (|45| . we 
obtain 



Finally the hopping parameters of both models are 
given in the Table IIIII 

Remark: 

In the WE model, we use E m as model parameter, 
hence T w e, and in the Harrison model, we use e m as 



a 


V s da 


Vpdcr 


V a 




3.9 


2.9855 


2.7237 


2.0754 


1.5590 


4.0 


2.7054 


2.4064 


1.8486 


1.3854 



Table HI: parameters of WE and Harrison models 




Ed + 6 



contribution of the CF acting on the e.g. p-states at the 
T point: The interactions with Sr states, with core states, 
with Madelung potentials etc. Therefore, e p is rather a 
model parameter than the true atomic energy, E p , of a 
2p state. If we speak about the model only, we may drop 
E p and Et2g, and retain only e p , Ed and E eg . 
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The parameters and V a are from the WE (p-d) model 
and Vpdvr — Vtt and V v d<y are from the Harrison (s-p-d) 
model, see App. I VI CI 

e denotes the energy of the atomic level and E, as used 
before, denotes the energy level corrected by the crystal 
field: A CT = E d - E p = e d - e p + 8 C f,<t, cf. Eq. (gSJ. Note, 



that ScF,m is different for m = n and m = cr, since Ed is 
the atomic energy level, and thus does not depend on m. 
This is the main reason, that, 5cf,o > &cf,k- Furthermore, 
ScF,a has strong dependence on a. 

Note, that we use another sign in the definition of n Ew 
compared to^ Eq (46) and (47). 



